function [mesh,row_idx] = smooth_C1_dis(dofs,mesh,cr)
% function [H,row_idx] = smooth_C1_dis(dofs,V,T,E,posE,ET,TE,d,cr)
% get the global smoothness conditionsidx_mat = zeros(3,d,3); % the column idx 1 2 3 -1 -2 -3
d=mesh.d;
idx_mat = zeros(3,d,3);
idx_c0 = zeros(d+1,3);
V=mesh.V;
T=mesh.T;
E=mesh.E;
posE=size(E,1);
ET=mesh.ET;
TE=mesh.TE;
d=mesh.d;
for i = 1:3
    idx_c0(:,i) = cr_indices(0,d,i,cr);
    idx_mat(i,:,1) = cr_indices(1,d,i,cr);
    line = cr_indices(0,d,i,cr);
    idx_mat(i,:,2) = line(1:d);
    idx_mat(i,:,3) = line(2:d+1);
end
E_in = find(ET(:,2)~=0); nE_in = length(E_in);
row_idx = zeros(nE_in*(2*d+1),1);
i = zeros(nE_in*(4*d+1),1); j = zeros(nE_in*(4*d+1),1); s = zeros(nE_in*(4*d+1),1);
pos = 1; row = 1;
v_idx = [1 2 3;2 1 3;3 1 2];
for k = 1:nE_in
    eg = E_in(k);
    tri = ET(eg,1); nei = ET(eg,2); % find two triangles who have eg.
    eg2tri = find(TE(tri,:)==eg);  % local index of eg in tri
    eg2nei = find(TE(nei,:)==eg); % local index of nei in tri
    % first set c0 condition:
    i(pos:pos+d) = row:row+d;
    j(pos:pos+d) = dofs(tri,idx_c0(:,eg2tri))';
    s(pos:pos+d) = 1; pos = pos + d + 1;
    % second set c1 condition:comput the barycentric of v4 relative to tri
    v1 = v_idx(eg2tri,1);    v2 = v_idx(eg2tri,2);    v3 = v_idx(eg2tri,3);
    A = [1 1 1;V(T(tri,v1),1) V(T(tri,v2),1) V(T(tri,v3),1);V(T(tri,v1),2) V(T(tri,v2),2) V(T(tri,v3),2)];
    b = A\[1;V(T(nei,eg2nei),1);V(T(nei,eg2nei),2)];
    % for each coordinates set c1 condition
    for l = 1:3 % the idx is always from v2->v3
        i(pos:pos+d-1) = row+d+1:row+2*d;
        j(pos:pos+d-1) = dofs(tri,idx_mat(eg2tri,:,l))';
        s(pos:pos+d-1) = b(l);
        pos = pos + d;
    end
    if mod(eg2tri+eg2nei,2)==1 % the same direction with v2->v3
        row_idx(row:row+d) = dofs(nei,idx_c0(:,eg2nei))';
        row_idx(row+d+1:row+2*d) = dofs(nei,idx_mat(eg2nei,:,1))';
    else  % the different direction in neighbour
        row_idx(row+d:-1:row) = dofs(nei,idx_c0(:,eg2nei))';
        row_idx(row+2*d:-1:row+d+1) = dofs(nei,idx_mat(eg2nei,:,1))';
    end
    row = row + 2*d + 1; 
end
dim = max(max(dofs))-min(min(dofs))+1; pos = pos-1; row = row-1;
H = sparse(i(1:pos),j(1:pos),s(1:pos),row,dim)-sparse((1:row)',row_idx,ones(row,1),row,dim); 
mesh.H=H;